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Abstract The efficacy of fast/slow MHD mode conversion in the surface layers of 
sunspots has been demonstrated over recent years using a number of modelling tech- 
niques, including ray theory, perturbation theory, differential eigensystem analysis, 
and direct numerical simulation. These show that significant energy may be trans- 
ferred between the fast and slow modes in the neighbourhood of the equipartition 
layer where the Alfven and sound speeds coincide. However, most of the models so far 
have been two dimensional. In three dimensions the Alfven wave may couple to the 
magnetoacoustic waves with important implications for energy loss from helioseismic 
modes and for oscillations in the atmosphere above the spot. In this paper, we carry 
out a numerical "scattering experiment", placing an acoustic driver 4 Mm below 
the solar surface and monitoring the acoustic and Alfvenic wave energy flux high 
in an isothermal atmosphere placed above it. These calculations indeed show that 
energy conversion to upward travelling Alfven waves can be substantial, in many 
cases exceeding loss to slow (acoustic) waves. Typically, at penumbral magnetic field 
strengths, the strongest Alfven fluxes are produced when the field is inclined 30°- 
40° from the vertical, with the vertical plane of wave propagation offset from the 
vertical plane containing field lines by some 60°- 80°. 
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1. Introduction 

Schunker and Cally (2006) (hereafter SC) and Cally (2007) recently presented a 
generalization of magnetohydrodynamic (MHD) ray theory that describes mode 
transmission and conversion between fast and slow magnetoacoustic waves in the 
context of solar active regions. A uniform, inclined magnetic field was imposed on a 
plane parallel atmosphere, and magnetoacoustic rays in the vertical plane containing 
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the field lines were modelled using the general ray transmission/conversion formalism 
of Tracy, Kaufman, and Brizard (2003). The ray results were fully confirmed by 
numerical solution of the wave equations. By remaining in the plane defined by the 
magnetic field and gravity though, only the magnetoacoustic modes are involved. 
The Alfven wave, the plasma displacement of which is perpendicular to this plane, 
is entirely decoupled. The purpose of this paper is to extend the wave analysis to 
the full three dimensions (3D), thereby bringing the Alfven wave into play. This 
occurs because the Alfven wave displacement vector is rotated out of the horizontal 
plane, and therefore interacts with the atmosphere's stratification. The situation 
differs from the coronal analysis of Melrose (1977) and Melrose and Simpson (1977), 
where field-line curvature and twist are responsible for tying the magnetoacoustic 
and Alfven waves together. 

Importantly, we integrate from interior to high atmosphere, so as to explore the 
coupling of helioseismic waves (p-modes) with atmospheric oscillations, especially 
in the vicinity of the Alfven/sound speed equipartition level z eq (typically shallow 
subsurface in sunspots) where fast/slow conversion occurs. 

Several important lessons may be drawn from our analysis. 

1. The 2D result (SC) that the atmospheric slow wave flux has a pronounced max- 
imum at a particular value of magnetic field inclination from the vertical 9, 
(typically in the range 20°-30° depending on spherical harmonic degree I), sur- 
vives into 3D, at least out till \<f>\< 90°, where <f> is the angle between the vertical 
plane of wave propagation and the vertical plane in which the field lines lie. 

2. There is a distinct maximum in upward Alfvenic energy flux in the atmosphere at 
around 9 = 40°, <f> = 60°, with these numbers only weakly dependent on model 
parameters. The Alfven wave can only be produced by mode conversion in this 
model. 

3. Conversion of magnetoacoustic waves to Alfven waves happens over an extended 
height range, typically some hundreds of kilometres, as distinct from fast/slow 
conversion which is normally localized sharply near z eq . 

In the interests of brevity, the reader is referred to SC for an extensive discussion 
of the solar context and significance to local helioseismology of MHD waves in strong 
magnetic fields near the surface, and to several earlier papers exploring MHD mode 
conversion using a variety of techniques (Spruit and Bogdan, 1992; Cally and Bogdan, 
1993: Cally, Bogdan, and Zweibel, 1994; Bogdan and Cally, 1997; Cally and Bogdan, 
1997;Cally, 2000; Crouch and Cally, 2003, 2005). 

A 3D analysis from the ray perspective is deferred until Paper II. 

2. Wave Mechanical Formulation 

The linearized wave equations and the necessary boundary conditions are introduced 
in this section. 

2.1. Equations 

In Section 3, a broadly realistic solar interior model will be coupled to an overlying 
isothermal atmosphere. A uniform, inclined magnetic field 

Bq = Bo (sin 9 cos <j), sin 9 sin <fi, cos#) , (1) 
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will be assumed throughout, where 9 < 90° is the field inclination from the vertical, 
and (f> the angle by which it is rotated from the x-z plane. The isothermal top 
adequately represents the chromosphere for our purposes, and conveniently allows 
us to take advantage of known exact series solutions when selecting top boundary 
conditions (see Section 2.2 and Appendix A). A horizontal and time dependence 
exp[i(fcx — cot)] is assumed, with z dependence to be determined, corresponding to 
wave propagation in the x-z plane. 

The linearized, adiabatic MHD equations for this scenario may be expressed in 
terms of the components of the displacement vector £ = (f , 77, £), the sound speed 
c, Alfven speed a, gravitational acceleration g, and density scale height H: 



2 2/2 2 2 \ ^ 

cos (j) sin 9 sin <fi r\ k — I cos 9 + sin 9 sin <j> J £ k 
+ sin 9 (i k sin (sin 9 sin 4>C' ~ cos $ v') ~ cos $ cos 4> C") + cos2 $ £" 

+ k ( a 2 k cos 9 cos 4> sin 9 -ig)C+ (uj 2 - ck 2 ) £ + ickC,' = 0, (2) 



i + — ifc cos sin 2 9 sin 77' 'a 2 +ik sin 2 9 sin 2 £ 'a 2 + (a 2 sin 2 f? + c 2 )£" 

ff 2 ^ / 2 2 2 2 2 \ 

- cos 9 sin f? sin <j>rj a — cos 9 cos sin 9 £ a + Iw — a k cos sin" 9 J £, 

fc cos # cos sin f?a 2 + i + (c 2 )'^ + (c 2 )'^' + i c 2 k £' = 0, (3) 



+ k£ 



and 



/2 Il2i2 \ 2 

i k cos f? £ c + cos £ c + I ai cos 9 — igk cos sin ) £ + w sin 9 sin 77 



i k cos 9 c . 2 2, 2x / • /j 1 • 1 a 2 

+ 1 <7K cos 9 + {u> — c k ) cos sm 9 + lk cos 9 c 



II 



, . , 2 cos 9 c 

+ [ 1 k cos sin 9 c — h cos 

H 



9c 2 '^(' = 0. (4) 



Primes indicate derivatives with respect to z. 

2.2. Fast, Slow, and Alfven Waves and Radiation Boundary Conditions 

The top boundary conditions for the fast and slow waves in the superposed isothermal 
atmosphere are straightforward. The asymptotic controlling factors for the displace- 
ments as z — > 00 are respectively (see the Appendix) 



exp[±fcz] and exp 



1 — 2ik H t&n9 cos ± i\ / ^ sec 2 9 — 1 ] 

4 ) 2H 



(5) 



where lo c = c/2H is the acoustic cutoff frequency. 

The fast wave is clearly evanescent, with exp[— kz] the appropriate solution. The 
term 2ikH tan 9 cos in the slow wave controlling factor is purely geometric, ac- 
counting for the change in x along an inclined field line: x — z tan 9 cos = constant 
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along a field line so the exp[ifcx] dependence cancels this term (recall that the slow 
wave is rigidly channelled along B in the high-altitude limit). The exp[z/2iJ] term 
is extinguished by the p 1 ^ 2 factors in the kinetic energy density ^Ip 1 ^ 2 ^ 2 , and 
similarly in the acoustic energy. The square-root term however is more important. 
If w > 0J C cos 0, the slow wave can propagate, and the "+" sign is the appropriate 
choice for it to be upgoing at infinity. On the other hand, for u> < u c cos6, i.e., 
below the ramp-modified acoustic cutoff frequency, the slow wave is evanescent. In 
that case too, the "+" sign is on the correct choice. 

The selection of appropriate top boundary condition for Alfven waves is more 
subtle. The difficulty is best illustrated by reference to the well-known exact solution 
for the perpendicular (to both B and g) velocity 77 in terms of Bessel functions in 
the 2D isothermal case <f> = 0: 

77 = s 2iKtane [A Jo (2s sec 6) + BY {2ssec9)} , (6) 

where s = uiH/a = (ujH/ao)exp(—z/2H), k = kH is a dimensionless horizontal 
wavenumber, and A and B are integration constants. If the atmosphere is deemed 
to extend to z = 00 (i.e., s = 0), the common practice is to set B = 0, resulting 
in a standing wave, the Jo solution alone (Ferraro and Plumpton, 1958; An et at, 
1989). With this choice, the Alfven wave may not take away energy, contrary to 
expectations. However, although the Yq Bessel function diverges as s — > + (z — > 00), 
its energy density is bounded, and so there is no need to dispense with it. The correct 
choice is actually B = —iA, whence 

v = As 2iKtane i4 2) (2ssec6>), (7) 
(2) 

(Schwartz, Cally, and Bel, 1984), where Hq = Jo — iYq is the second Hankel 

function of order zero, which is well known to represent a wave travelling in the 

negative s (positive z) direction. Although the Alfven wave (travelling at speed a = 

ao exp(z/2H)) actually reaches z = 00 in a finite time (2H/ao from z = 0), it is what 

it does there that is important: if it is assumed that it reflects totally, a standing wave 

quickly results (the Jo solution), but if on the other hand it is assumed lost there, the 
(2) 

solution develops. Of course, in reality the isothermal atmosphere (representing 
the chromosphere) does not extend to infinity. Nevertheless, if we imagine that the 
wave is lost once it reaches a great height above the solar surface, which is the 
more physical scenario, the Hankel solution is appropriate. Because we are focussed 
on the conversion process per se, we choose not to complicate issues by including 
a multi-layer (chromosphere/corona) atmosphere and the resulting resonant partial 
reflection from the transition region (Cally, 1983; Schwartz, Cally, and Bel, 1984). 

The full three-dimensional case is more complicated, because it inherently involves 
coupling between all three MHD modes: Alfven, fast, and slow. Closed form solutions 
in terms of special functions are not available. Instead, a Frobenius series expansion is 
required. This is carried out in the Appendix, where the construction of the outward 
and inward travelling wave solutions is explained. 

With all of this in place, the three top boundary conditions in the numerical 
solution of the full 3D sixth-order coupled wave equations consist of a matching onto 
the three physical solutions in the isothermal atmosphere: the evanescent fast wave, 
the outgoing or evanescent slow wave, and the outgoing Alfven wave. 
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Figure 1. Acoustic (full curve) and magnetic (dashed) wave energy fluxes as functions of 
height (z) over the full computational domain, for two different frequencies: 3 mHz (left panels) 
and 5 mHz (right panels). In all cases, the 2 kG magnetic field is inclined at 9 = 30° to the 
vertical. The top panels refer to the 2D case (<f> = 0°), the middle panels to <j> = 30°, and 
the bottom panels to <j> = 60°. The acoustic cavity base is at z\ = —5 Mm throughout, 
corresponding to k = 1.37 Mm -1 , or i fa 955. The vertical axis at z = —0.297 Mm indicates 
the position of the a = c equipartition. 



The bottom magnetic boundary conditions are applied deep enough in the interior 
that a <C c, and so both the slow and Alfven waves are rigidly field-guided and 
highly oscillatory with respect to the density scale height. This suggests the radiation 
condition 

(|—v)vx* = 0, (8) 

where a = a Bo is the Alfven velocity. The term in brackets is the field-directed 
upgoing wave operator, and it is applied to V x £ so as to suppress the acoustic 
wave (which is irrotational to a high degree of accuracy) and accentuate the highly 
oscillatory magnetic waves. 
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3. Numerical Model and Solutions 

We perform a "scattering experiment" (in the spirit of Spruit and Bogdan, 1992) by 
placing a plane acoustic source at Zb < beneath the solar surface. The Magnetically 
Modified Model S (MMMS), a variant on the well-known non-magnetic Model S of 
Christensen-Dalsgaard et al. (1996), is used below z = 0.5 Mm, and an isothermal 
slab above. This model, which reduces the gas pressure (and consequently the sound 
speed) near the surface to compensate for the introduction of a magnetic pressure, is 
detailed in the appendix of Schunker and Cally (2006). In practice, we use Zb = —4 
Mm, though this is not important, providing it is well below the acoustic - Alfvenie 
equipartition level z eq « —297 km (for 2 kG held) where a = c. At Zb, an acoustic 
wave of frequency uj is generated which has a natural cavity depth z\ < Zb- This 
determines its horizontal wavenumber (fc) or alternately spherical harmonic degree 
(£). Since a < cat zj,, the slow mode is predominantly magnetic and transverse, much 
like the Alfven wave. Radiation conditions are applied on both: no slow or Alfven 
waves enter the computational region from below. Together with the amplitude of the 
acoustic driver (set by normalizing the total acoustic energy per unit horizontal area 
in Zb < z < to unity), these constitute the three lower boundary conditions. The 
upper radiation - evanescence conditions are as prescribed in Section 2.2. Subject to 
these six boundary conditions (one of which is non-homogeneous), Equations (2-4) 
are solved numerically in Zb < z < zt- We are then free to adjust oj, z±, 9, <fi, Bq, 
etc. and to observe the vertical components of acoustic and magnetic wave energy 
fluxes F ac and F mag reaching the computational top zt, typically set at 2 Mm. 
For reference, we note that the vector wave energy flux is 

F = F ac + F mag = Re[piv*+exb*], (9) 

(Bray and Loughhead, 1974), where pi = — pc 2 V-£ + pgC 1S t fle Eulerian gas 
pressure perturbation, v = — is the plasma velocity, b = Vx (£xBo) is the 
magnetic field perturbation, and e = — vxBo is the electric field perturbation. 
The energy density is made up respectively of kinetic, acoustic, gravitational, and 
magnetic parts, 

E= P W f + \Plf + E N 2 K f (io) 

2 1 1 2pc 2 2 M 2/x ' y ' 

where N is the Brunt- Vaisala frequency. E and F satisfy the expected conservation 
equation dE/dt + V - F = by construction. 

Figure 1 displays the vertical fluxes as functions of z for various cases with a 2 
kG magnetic field inclined at 30° to the vertical. The top-left panel corresponds to 
a frequency of 3 mHz with cf> = 0°, so the Alfven wave is entirely decoupled, and 
cannot be excited by the acoustic driver at Zb = —4 Mm. All magnetic flux in this 
case is therefore associated with the magnetoacoustic waves, predominantly the slow 
wave below z eq , where it is negative, in accord with the lower radiation boundary 
condition. Since the (magnetically dominated) fast wave in z ;§> z eq is evanescent in 
all cases, as is the (acoustic) slow wave here since u> < to c cos 9 (the acoustic cutoff 
frequency is 5.2 mHz in the isothermal slab), both magnetic and acoustic fluxes 
vanish quickly as z increases. In the top-right panel, the frequency is increased to 5 
mHz, with other parameters unaltered. In this case, the ramp effect is sufficient to 
reduce ui c cos 9 below the wave frequency, and so the slow (acoustic) wave propagates 
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Figure 2. Acoustic (top) and magnetic (bottom) fluxes at z = zt = 2 Mm as a function of 
magnetic field angles 8 and </> for the 5 mHz, 2 kG case of Figure 1. Contours and shading 
indicate wave energy flux on a linear scale, with white being highest. 



upward indefinitely carrying an asymptotically constant positive flux. In the absence 
of the Alfven wave, the magnetic flux once again vanishes with increasing height. 

These two cases are repeated in the middle panels, but for <f> = 30°. Introduction 
of 3D coupling has brought in the Alfven wave, as is apparent from the constant 
positive magnetic flux high in the atmosphere. In the 3 mHz case, all upward flux 
is magnetic, since the acoustic wave is evanescent. At 5 mHz, both acoustic and 
Alfvenic flux contribute substantially. 

Finally, the bottom panels repeat the simulations with <fi = 60°, where the Alfven 
coupling is even stronger. 

Figure 2 illustrates the net outward acoustic and magnetic fluxes for 5 mHz waves 
in a 2 kG magnetic field of various inclinations (9) and angles of wave incidence (</>). 
As shown in SC and Cally (2007) in the case 4> = 0°, there is a prominent maximum 
in acoustic flux at around 9 = 26° (for z\ = —5 Mm). We now see that this extends 
out beyond 90° in <fi, though past say 50°, the height of the maximum reduces 
significantly. The surprising breadth of the maximum in the </> direction is discussed 
further in Section 4. Of course, there is no acoustic flux for cos 9 > uj/uj c (9 < 10.7°), 
because of the acoustic cutoff. 

On the other hand, Alfvenic (magnetic) power peaks at around 9 = 40°, <f> = 60°, 
and has a slightly larger maximum than does the acoustic flux. Very similar results 
are obtained for 1 kG (z eq = +53 km) and 3 kG (z eq = —613 km) magnetic field. 

Figure 3 again concerns 2 kG field, but for a deeper acoustic cavity (lower t), and 
for a lower frequency, with similar results. In the left panel, it is the acoustic flux 
that has the slightly higher maximum. In the 3 mHz case (right panel) there is no 
acoustic flux for 9 < 53.9° because of the acoustic cutoff, and very little thereafter, 
though the magnetic flux is broadly unchanged. 

By focussing on the fluxes alone, we have lost sight of the atmospheric fast wave, 
which is evanescent. However, we may redress this by examining the energy densities. 
Figure 4 displays these for a range of 9 and </> in the 5 mHz case of Figures 1 and 2. 
The first three panels display 2D cases (</> = 0°), where the Alfven wave is absent. 
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Figure 3. Same as Figure 2, but for different z\ (left) and frequency (right). Left: 5 mHz wave 
with acoustic cavity base z\ = —15 Mm, corresponding to k PS 0.67 Mm" 1 (I = 468). Right: 3 
mHz frequency with z\ = —5 Mm (k = 0.79 Mm -1 , I = 552). The acoustic and magnetic flux 
contours are shaded separately in the right panel, since the acoustic flux is much less than the 
magnetic flux. 



The magnetic energy density E mag above z eq is therefore a good indication of the 
strength of the fast wave there. With 9 = 0° (both fast and slow waves evanescent), 
E ac and E mag are of comparable magnitude, though much smaller than E^ n and 
Egrav- The slow wave is apparently more gravity wave in nature than sound wave; 
hardly surprising as the Brunt- Vaisala frequency in the isothermal slab is 4.9 mHz. 
At 9 = 30° the slow wave is travelling, and E ac 3> E mag in the atmosphere. In the 
highly inclined field case 9 = 60° though, where the theory of SC predicts that the 
coupling to the slow wave is weak, the magnetic energy density dominates between 
roughly 200 km and 800 km, despite fast wave evanescence. Finally, the 3D case 
9 = 40°, <p = 60°, corresponding approximately to the magnetic flux maximum, 
illustrates how the travelling Alfven wave can dominate the atmospheric energy 
density. Overall, these figures give a good indication of how much "leakage" of p-mode 
energy into the overlying atmosphere we can expect in active regions. 



4. Discussion and Conclusions 

The greater breadth of the acoustic maxima in the cf> direction compared to the 9 
direction exhibited in Figures 2 and 3 may be broadly explained using the concept 
of attack angle discussed in SC, i.e., the angle a between the wavevector k and 
the magnetic field Bo- Ignoring the acoustic cutoff and Brunt- Vaisala frequencies, 
generalized ray theory predicts an acoustic transmission coefficient of 



exp 



-7T |k| h sec ijj sin 2 a 



(11) 



where tp is the angle the wavevector makes to the vertical at the equipartition level 
z eq , and h = [d{a 2 / c 2 ) / dz\a= c is the thickness of the a fa c layer. In 3D, sin 2 a = 
1 — (sinacos^sin-i/; + cos^cos';/') 2 - If we crudely set |k| = co/c, then for the case 
of Figure 2, |k| = 3.7 Mm -1 and h = 0.29 Mm, i.e., kh w 1.1. Figure 5 presents 
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Figure 4. Acoustic (full curve), kinetic (dotted), magnetic (short dashed), and gravitational 
(long dashed) energy densities for the 5 mHz case of Figure 2, with 8 and <j> as labelled. 



contours of T in <$> — Q space for this case, with ip = 22°, and clearly displays a greater 
breadth in <f> than in 9. If the acoustic cutoff effect is reintroduced, the corresponding 
flux contours must be squeezed upwards at the bottom to produce a "blank zone" 
below 9 = 10.7° (as in Figure 2) and an even more pronounced anisotropy. This is 
very much in qualitative agreement with the acoustic results of Figures 2 and 3. 

It is clear from the few examples presented in Section 3 that 3D magnetoacoustic 
coupling to the Alfven wave is very strong, rivalling and in some cases surpassing the 
maximum fast-to-slow transmission. Once the angle (</>) between the vertical plane 
of the wave and the vertical plane of the magnetic field differs, the coupling turns on. 
It is typically strongest at <f> = 50°- 70°, and for 9 = 30°- 40°. This is a significant 
result, with important implications for our understanding of p-mode lifetimes, and 
atmospheric oscillations above sunspots. 

The assumption of uniform magnetic field is made for simplicity and in the spirit 
of a simple scattering experiment. It is probably quite reasonable for the localized 
fast -slow mode conversion. However, the greater breadth of the magnetoacoustic - 
Alfven coupling region suggests that field-line curvature may be relevant to that 
process. Curvature can certainly provide another Alfven coupling mechanism in its 
own right. 

The upward (Alfven) and downward (slow plus Alfven) magnetic losses illustrated 
in Figure 1 indicate a range of relative efficiencies. At fixed 9 and <f> though, it seems 
that upward loss through the atmosphere is more important at higher frequency. 
And of course, upward loss is typically of greater relative importance in that 9 — (f> 
region where F mag is maximal. 

The model presented here differs from the earlier Crouch and Cally (2005) analysis 
of the interaction of helioseismic p-modes with 3D inclined magnetic field in that (i) 
it consists of a realistic (albeit modified) solar interior model rather than a polytrope; 
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Figure 5. Transmission coefficient T in the <fi— 6 plane as given by Equation (11) with 
\k\h = 1.1 and ip = 22°. The contour levels are 0.9, 0.8, 0.7, .... 

(ii) it contains an overlying atmosphere, rather than a vacuum, which allows escape 
of waves upward as well as downward; and (m) it includes a driver instead of solving 
for (undriven) eigenmodes. Since our focus is on coupling to atmospheric oscillations, 
point (ii) is especially important. Nevertheless, it is interesting that that study also 
found substantial (downward propagating) Alfven wave coupling. 

Only the relative magnitudes of the acoustic and magnetic fluxes were quantified 
in Section 3. In practical terms though, how big are these fluxes, and how quickly 
can they deplete the underlying p-modes which drive them? Recall that the fluxes 
were normalized by setting £ = f E ac dz to 1, with z^ = —4 Mm in all cases. Now, 
F I £ has the dimensions of inverse time, and broadly represents the inverse time scale 
over which the p-mode energy in z^ < z < depletes, assuming it were not being 
continually replenished by the driver at Zb- The maximum fluxes seen in Figures 2 
and 3 are of order three in these normalized terms, where time is being measured 
in kiloseconds. Consequently, the minimum decay time is of order t ~ g ks « 300 
seconds only, i.e., not much more than one period! Of course, this is somewhat of an 
underestimate of r, because: 

1. it relies on a particular combination of 6 and </>; only waves travelling in a par- 
ticular direction and meeting field of a particular inclination will be so strongly 
affected; 

2. we have not integrated the acoustic energy over the entire acoustic cavity z\ < 
z < 0, so the available energy store is larger than indicated by £ ; 

3. there is also a substantial store of kinetic energy in the acoustic cavity on which 
to draw. 

Nevertheless, the surface magnetic field and its attendant mode conversion (to slow 
or Alfven waves) is clearly significant for p-modes in active regions. 

Figure 1 helps us understand the nature of the mode conversion. In all cases, 
we see a broad conversion region starting a little below z eq , but extending well 
above it. This is quite different from the fast -slow transmission - conversion process, 
which is strongly localized around the equipartition level. The reason for this may 
be discerned using a WKB description, which will be addressed in the forthcoming 
Paper II. 
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Appendix 

A. Frobenius Solutions in an Isothermal Atmosphere 

Exact solutions for the MHD wave equations in an isothermal atmosphere with 
uniform magnetic field were developed by Zhugzhda and Dzhalilov (1984). In the case 
where the wave vector lies in the vertical plane containing the field, the magnetoa- 
coustic waves decouple from the Alfven wave, and were expressed in terms of Meijer 
G-functions (although Cally, 2001, showed later that the simpler 2-F3 hypergeometric 
functions could be used instead). However, in the full 3D case, solution in terms of 
special functions does not appear to be possible, and Frobenius series were developed 
instead. 

When the two magnetoacoustic waves are coupled to the Alfven waves, the gov- 
erning differential equations are of sixth order, but there are only five distinct roots 
of the indicial equation. The Alfvenic solutions display a double root, and hence 
one logarithmic solution. We need this logarithmic solution in order to represent the 
outgoing Alfven wave. Unfortunately, Zhugzhda and Dzhalilov (1984) did not pursue 
the logarithmic case. We therefore develop the necessary solution here. 
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We adopt the same dimensionless variables as in the case of the decoupled Alfven 
wave, v = uH/c, s = toH/a = (ujH/aQ)exp(—z/2H), and k = kH. The MHD 
wave equations can be represented as a set of three coupled second-order ODEs 
in the displacements £(s), r](s), and C( s ); derived from Equations (2-4). However, 
computationally, we find it more convenient to write them as a single sixth-order 



matrix equation 

sU' = AU, 

where U(s) = C, a£', sr, 1 , sff = (£ , s £') T and 



(12) 



I 

P QJ ' 



(13) 



and where each of the four constituent blocks is 3 x 3. Specifically, 



Q 



2i« cos (f> tan ( 


2iK 



-2i« sin d> tan 9 2ire 



4i« cos (j> tan ( 




and the components of P are 



se f-Z + tan 2 I 



2 cos d> tan 9 



—2 sin (j> taxi 6 
2iKCOs<Atan# — 2 



(14) 



Pn 

Pl2 
Pl3 

P21 
P22 

P23 

P31 
P32 
P33 



tan 



1 / + 



4k' 



4 ) s z + 4k - 2z/ sec 



-2 (cos2</> 

+ 4i(l — 7 _1 )KCOs0tan 
— 2(k 2 + v 2 ) sin2</>tan 2 9 

4i (cos 4> tan 9 (17 (k 2 + v 2 ) + Kcos^tan^) i/ 2 + s 2 Ksec 2 9) 
4i sin (j) tan (17 cos 4> tan #j/ 2 + (7 — 1)k) 



7 

a ( 2 2 Q . ( 2 . 2 

4 I s sec & + I v sin 



k cos <f> J tan 



4 sin (j) tan (7V 2 — iK cos tan (?) 



4i (7-l)^ +4(K 2_^ )cos ^ tang 

7 

—Au 2 sin^>tan(? 
4i« cos«!)tan( 



4z/ 



7 



(15) 
(16) 

(17) 

(18) 
(19) 

(20) 

(21) 
(22) 
(23) 



where 7 is the ratio of specific heats. Note that A(s) = Aq + A2S 2 



A partial theory of Frobenius expansion for matrix equations is set out in Rubinstein (1969). 



ray3DIrev.tex; 2/02/2008; 7:12; p. 12 



3D MHD Wave Propagation I 



Adopting the standard Frobenius expansion about the regular singular point s = 
(z = oo), 

OO 

U(a) = U « sn+M . ( 24 ) 

n=0 

the lowest order balance yields the indieial equation 

A u = /m , (25) 

indicating that fi is an eigenvalue of Ao and uo is the corresponding eigenvector. 
The spectrum is 

(i g |— 2k, 2k, —iy/ 4i/ 2 — cos 2 6 sect? + 2iKCOs^tanfl — 1, 

\\J Av 2 — cos 2 6 seed + 2ik cos <j> tan 6 — 1, 2iKCOs</>tan6', 2iKCOs^tan6'|, (26) 

corresponding respectively to the growing (unphysical) fast mode, the evanescent 
fast mode, the outgoing slow mode (assuming v > \ cos0), the incoming slow mode, 
and the Alfven mode. As mentioned before, the Alfven eigenvalue has algebraic 
multiplicity 2, although its geometric multiplicity (dimensionality of its eigenspace) 
is only 1. 

Beyond the first order, we find that ui = in all cases, and derive the recurrence 
relation 

u„ = - (A - (n + v)\y x A2U„_ 2 , (27) 



from which it is apparent that all odd coefficients vanish. This completes the de- 
scription of the first five solutions. 

To discover the sixth solution, we must replace (24) with 



U 6 (s) 



U 5 (s) In s + v„s n 



+A»6 



(28) 



n=0 



The factor 2/n is just a convenient normalization. Substituting this into (12) and 
equating coefficients we find 



(A - M6l)v = u (29) 
vi = (30) 
(A - (n + /x 6 )I) v„ = u„ - A 2 v„_ 2 , (31) 



where the u„ are the coefficients in the first Alfvenic solution U5. Equation (29) 
yields a generalized eigenvector, in the sense that Ao — nel is singular by definition, 



ray3DIrev.tex; 2/02/2008; 7:12; p. 13 



P. S. Cally and M. Goossens 



with uo in its null space, i.e., (Ao — /i6l) 2 vo = 0. We find 

i cos 3 9 sin 9 sin <j) cos 2 9 i cos 2 9 cos cj) sin 2 9 sin < 



vo 



2k (l - sin 2 9 sin 2 4>) ' 2jis 2 ' 2k (sin 2 9 sin 2 - l) 



sin 2 6> ((cos 20 - 3) sin 2 + 4) sin 20 2 2 i K CO s 6> cos 4> sin 

: , 1 - sin 0sin <p H 5 > 

4 sir 6»suT - 1) 7^ 



, \ T 

cos sin 6 sin <f> \ 

sin 2 sin 2 6 - 1 / 



(32) 



Naturally, this is determined only up to an arbitrary multiple of the Alfvenic eigen- 
vector 



uo = ^— cos 0sin 2 9 sin0, 1— sin 2 9 sin 2 cf>, — cos # sin # sin 0, 

— sin 2 sin' 

2 \ 1 

2iKCOs^sin #sin0l , (33) 



2i«;cos 2 0sin 2 #sin0tan#, 2irecos0(l — sin 2 #sin 2 <fi) tan#, 



though (32) has been engineered to yield a solution U6 with zero net energy flux, 
i.e., it is a standing wave, as we shall see shortly. 

The remaining task is to determine which combinations of U5 and Uq rep- 
resent incoming and outgoing Alfven waves at z = 00. To do this, we evaluate 
the field-aligned component of the Poynting flux, Fn = Bo-F maff , where F ma g = 
Re[io;(^xBo) x Vx (£* x Bo)], as s — > + , which only involves the n = coefficients. 
Letting U = C5U5 + iCeU6, we find 

F|| = -|F cos 9 (CfC 6 + C 5 Ce) (l - sin 2 9 sin 2 <pj , (34) 

where Fq = B 2 ci//(nH 2 ). Notice that F\\ = if either C5 = or C6 = 0, indicating 
that U5 and U6 are both standing waves. Also note that F\\ is invariant under the 
transformation C5 — > C5 + i a Cq, aGl, and similarly for Cq. 

Setting C5 = A+ + A- and 6*6 = A_ — we may alternatively represent the 
general Alfven wave as a linear combination of upgoing and downcoming modes, 
U = A+XJ+ + A-XJ-, where U± = (1 =F i«)U5 =p iU6, with total flux 



*|| = F 



|A+| 2 -U_| 2 



cos Oil — sin 2 9 sin 2 <f> ) ■ (35) 



The vertical flux is F z = Fn cos 9. The real coefficient a does not affect the energy 
flux, but does contribute to the energy density and the solution matchings. To 
determine it, we look for guidance in the two dimensional case. 
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A.l. 2D Case 

Consider the 2D case = 0, where £ = £ = and only the transverse displacement 
77 remains: 



2iK tan 6 T / 

775 = s Jo(2ssec t 



(36) 



Y (2ssec9) + 



cos 9 2 In sec ( 



2C 

— I Jo (2s sec 9) 



(37) 



where C = 0.577216... is Euler's constant. Clearly, Jo(2ssec0) and Yb(2ssecf?) 
are linearly independent solutions, as expected. The corresponding outgoing wave 
solution then takes the familiar Hankel function form 



77+ = s 2iKtanfl [(1 _ i/3)J (2ssec<9) - i Y (2s sec 9)} 



2iK tan t 



H^ } (2ssec9) - i/3 J (2s sec ( 



(38) 
(39) 



where 3 = a + (^4 ~ 21nscce - — ). It is well-known 2 that the H {2) Hankel 
function represents a "pure" wave travelling in the negative s (positive z) direction, 
so it is apparent that we must set 8 = 0. This determines a in the 2D case. 



A.2. 3D Case 



In 3D, we do not have the luxury of closed- form solutions from which the large s 
asymptotic behaviour may be determined. In fact, this would not even be appropri- 
ate, as the Alfven wave undergoes coupling to the magnetoacoustic waves around 
s = 0(1). We are trying to find the real a which delivers a pure outgoing wave at 
infinity, despite it not exhibiting sinusoidal behaviour there. But how do we identify 
a "pure" wave in this regime? To do this, we can again use the Hankel function, or 
at least its asymptotic behaviour as s — > + , 

H^\2ssec9) ~ 1 - — (Ins + In sec 6? + C) + 0{s 2 Ins). (40) 

(As s sec 9 is distance along a field line, independent of 0, it is the appropriate spatial 
coordinate for the field-guided Alfven wave.) A pure outgoing wave should have this 
asymptotic structure. 

Now, it is easily confirmed that the asymptotic polarization direction of £ + as 
s — > is d = (— cos 0sin 2 9 sin0, 1 — sin 2 9 sin 2 0, — cos 9 sin 9 sin 0) , which is 
perpendicular to Bo as expected, and that the small s behaviour of displacement in 
this direction is 

, j. . . 2 i i cos 2 9 . . 

d-£, ~1_ i a Ins (41) 

7T TT"/l> Z 



2 This is verified using the large argument asymptotic behaviour (Abramowitz and Stcgun, 
1965), H ( 2 \x) ~ y / 2/(wx)exp[-i(x - tt/4)]. 
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where an arbitrary normalization has been suppressed. Comparing this with Equa- 
tion (40), we again infer that 

2 In sec 9 2C cos 2 6 , t . 

a = + j, 42) 

as in 2D. 

In summary, the three required radiation solutions at large z are U2, U3, and 
U+ = (1 — ia)U5 — iUg. Physical solutions must match to a linear combination 
of these. In practice, we apply the boundary conditions at z = 2 Mm, where s = 
uH/a ~ 0(1O~ 4 ) or smaller, so convergence of the series is very rapid. 
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